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Abstract In this paper, a class of deterministic sensing matrices are con¬ 
structed by selecting rows from Fourier matrices. These matrices have better 
performance in sparse recovery than random partial Fourier matrices. The 
coherence and restricted isometry property of these matrices are given to eval¬ 
uate their capacity as compressive sensing matrices. In general, compressed 
sensing requires random sampling in data acquisition, which is difficult to im¬ 
plement in hardware. By using these sensing matrices in harmonic detection, a 
deterministic sampling method is provided. The frequencies and amplitudes of 
the harmonic components are estimated from under-sampled data. The sim¬ 
ulations show that this under-sampled method is feasible and valid in noisy 
environments. 

Keywords Compressed sensing • Deterministic sensing matrices • Under¬ 
sampling • Harmonic detection 


1 Introduction 

Compressed sensing (CS) theory asserts that one can recover sparse signals 
from far fewer samples or measurements than traditional methods [3. The 
concept of CS is to recover high dimension sparse signals through low dimen¬ 
sion measurements. The two fundamental questions in compressed sensing are 
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how to construct suitable sensing matrices and how to recover sparse signals 
from under-sampled data efficiently. To find the sparsest solution, one would 
solve the problem 

argmin||x||o s.t. 4>x = y, ( 1 ) 

where x G is a sparse vector, Hq norm denotes its number of nonzero 
elements, $ G is called a (compressed) sensing matrix and y G 

is the measurement. However, this is a hard combinatorial problem. In order 
to make problem © solvable, the sensing matrix $ must obey a uniform 
uncertainty principle. 

The sensing matrix $ has the {k, S)-restricted isometry property (RIP) if 

(l-(5)||x|l^<||$x||^<(l + <5)||x||^ (2) 

holds for all fc-sparse vectors x, ||x ||2 denotes .^ 2 -norm of x [5]. The smallest 
<5 for (fc, 5)-RIP is the restricted isometry constant (RIC) Sk- If a complex 
matrix $ satisfies (fc, 5)-RIP, then the 2M x 2N real matrix formed by 

replacing each element a + \/—lb by matrix ^also satisfies the same 

(fc, ^)-RIP in real number domain [3]. Let # be a matrix with £ 2 -normalized 
columns j^Pn, i-e. ||y >„||2 = 1 for n = 1,2, ■ ■ ■ , N, the condition 

© is equivalent to that the Gram matrix of every column submatrix 

4>k:(/C C {1,2,-- - ,A^},|/C| < k) has all its eigenvalues in the interval [1 — 
5k, 1 + ^fe]. If 52k < 1, the problem ([T]) has an unique Ic-sparse solution. It has 
been proven that when 52k < — 1 the problem ([T]) can be approximated by 

a relaxed £i norm convex optimization [4] , that is 

argmin||x|| s.t. $x = y. (3) 

Another criterion often used to evaluate the property of a CS matrix is 
coherence. For the matrix €> with ^ 2 -normalized columns, the coherence of $ 
is defined as 

= max \{pn,Pn')\- (4) 

l<n^n'<N 

Although the computation of coherence which only involves two columns each 
time is much more feasible, it considers the worst case and often leads to 
results somewhat pessimistic. 

By the concentration inequalities, some random matrices are proven to 
satisfy RIP with high probability, such as matrices with independent Gaussian 
or Bernoulli elements [T], or matrices whose rows are randomly selected from 
the discrete Fourier transform matrices |12] . In practical application, random 
sensing matrices usually imply random sampling in data acquisition, which 
is difficult to implement in hardware. Many scholars have become interested 
in finding deterministic RIP matrix constructions mm- A deterministic 
construction of sensing matrices was given by polynomials over finite fields 
and the RIP weaker than random matrices was proven [7] . The class of partial 
Fourier matrices is of special importance in CS. Some deterministic sensing 
matrices of partial Fourier were constructed in and m- Sensing matrices 
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similar to partial Fourier matrices were also constructed explicitly m , even 
some matrices break the barrier related to RIP [3]. 

In this paper, we shall verify that a class of deterministic partial Fourier 
matrices can be used as sensing matrices and utilize these matrices to design a 
deterministic sampling method for harmonic detection. It will be demonstrated 
that the frequencies and amplitudes of harmonic components can be estimated 
from deterministic sub-Nyquist samples. 


2 Deterministic Sensing Matrices 

In this section, we give the sensing matrices directly and analyze their property 
as sensing matrices. Let Fa? be the N x N Fourier matrix whose element in 
TO-th row and n-th column is given by 

I jQi7r7nTL\ 

FN[m,n] =exp (— — — j , in,n = 1,2, ■■■, N, (5) 

where j = \/—l is imaginary unit. We restrict IV = 4z + 3 (z is a positive 
integer) to be a prime number and choose M = 2z + 1 = {N — I)/2 rows from 
Fjv to construct a.n M x N partial Fourier matrix. The indexes of the rows 
are given by 


M = {to ■. p ■ mod iV, to = I, 2, • • • , M} , (6 ) 

where p is an arbitrary positive integer co-prime to N . The M row indexes 
just don’t repeat according to number theory. Because of the periodicity of 
trigonometric functions, the elements of generated matrix A can be written 
as 

A[TO,n] =exp(^:^^^^|P^^ , m&M,n=l,2,--- ,N. (7) 


2.1 The Coherence 

Then we check the property of the matrix A. In order to calculate the coher¬ 
ence of A, we introduce the following theorem about quadratic Gauss sum: 


Theorem 1 (Theorem 1.5.2 in [2]) Let k and N be coprime integers with 
N > 0 and N = 3(mod4). T/ien 


N-l 


' (k; N)=Y^ exp i 

m—0 


/ j2'nkm? 


N 



where denotes the Jacobi symbol. 
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The coherence of A can be expressed as 

M 


a (A) = max — 


^expf 

m—1 


f j2Trdpm^ 


N 


( 8 ) 


Taking note that the values of {dprn^ mod N) distribute symmetrically when 
m = 1, 2, • • • , TV — 1, so we obtain 


p{A) = 


±jy/N-l 


2M 


VM+ 1 

V2M 


(9) 


It has been proven that orthogonal matching pursuit (OMP) and ba¬ 
sis pursuit (BP) both can recover the A:-sparse solution of ([T|) when k < 
So A can guarantee the recovery of sparse signals below the 

sparse level of ^ r; As what metioned in Sect.[Tl the result 

is somewhat pessimistic compared with numerical experiments. 


2.2 The RIP 


We give the RIP of A according to the following result: 

Theorem 2 (Theorem 14.1 in [5]) Suppose X is a k-sparse vector which is 
uniformly distributed among all k-sparse vectors. If for S,e G (0,1), there 
exists a constant c > 0 such that 


h{^)< 


cS 


In (A/e)’ A 


1 $ 


l2< 


cS'^ 

In (A/e) 


( 10 ) 


Then $ with £ 2 -normalized columns satisfies (k, 6)-RIP with probability at least 
1 — e. 


To present the statistical RIP of A intuitively, we plot the maximum and 
minimum eigenvalues of Gram matrices where are randomly se¬ 

lected column submatrices. We set M = 11 and A = 23. The Maximum and 
minimum eigenvalues of sub-Gram matrices ^ic are equivalent to 1 -I- (5 and 
1 — (5. These eigenvalues of a random partial Fourier matrix are also plotted 
for comparison. The data are obtained from 2000 sub-Gram matrices for each 
k. The solid lines sketch the average values of maximum and minimum eigen¬ 
values of all sub-Gram matrices and dashed lines sketch the limiting values. 
Fig. [1] shows that the eigenvalues of A’s sub-Gram matrices distribute more 
intensively than the random partial Fourier matrix. In fact, the correlation of 
any two columns of A is either {j'/N — l)/2 or — l)/2 with equal 

probability, which makes the variants of sub-Gram matrices far fewer. 

Then we plot the success probabilities of recovering fc-sparse vectors x via 
A and random partial Fourier matrices when M = 51 and A = 103. For 
convenience we use OMP algorithm. At every time, the locations of nonzero 
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(a) 



k 


(b) 

Fig. 1 Maximum and minimum eigenvalues of sub-Gram matrices for k. (a) A; (b) The 
random partial Fourier matrix. 


elements in vector x distribute uniformly among all fc-sparse vectors and the 
values of nonzero elements follow standard normal distribution. If the locations 
of nonzero elements in x can be found by OMP, we say that this trial succeeds. 
There are 10000 trials for every k. Fig.[2]shows that A has better performance 
than random partial Fourier matrices. 


3 The Application in Harmonic Detection 

Harmonic detection is widely used in electrical power systems [5]|11|[TB]. In 
m, a staggered under-sampling algorithm is proposed to detect the harmonic. 
The algorithm reorders sub-Nyquist samples before applying the fast Fourier 


















































6 


Shan Huang et al. 



Fig. 2 Comparison of recovery performance between A and random partial Fourier ma¬ 
trices. 


transform (FFT). We shall demonstrate that even fewer samples are enough by 
using CS. In this section, we analyze the problem and derive our under-sampled 
method in complex number domain. Then the practical implementation in 
real number domain is given later. The following discussion is in noiseless 
situations, but the method is also robust to noise as simulations show. 


3.1 Problem Formulation 

The signal s{t) contains K harmonic frequency components and we intend to 
estimate the frequencies and corresponding amplitudes. The samples can be 
expressed as the sums of K complex exponentials, namely 


K 

s{t) = ^ Cfc exp (j (27r/fc • t + Ok)), 

fc=i 


( 11 ) 


where Ck, fk and 6k represent the amplitude, frequency and phase angle of 
the fc-th harmonic component, respectively. The frequencies are all integer 
multiples of some fundamental frequency /o, i.e. fk € {/o, 2/o, • • • , N fo}. 

Suppose we have an analog-digital converter with sampling rate /s = 1 /At, 
arbitrary samples at time points At, 2At, ■ ■ ■ , I At, ■ ■ ■ can be obtained. In 
general, the sampling rate fs must be higher than the Nyquist rate. But 
because the harmonic signal is sparse in frequency domain, we can use a specific 
low-rate analog-to-digital converter to detect harmonic components. 
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3.2 The Method for Under-sampling 

Since the harmonic components of the signal are aligned with a discrete fre¬ 
quency grid, the sample s[l] at time point lAt can be written as 

N 

S W = ^ Cn exp (j (2TTfn ■ I At ^ (12) 

n—1 

Certainly, without regard to noise, Cn corresponding to /„ = n/o (n = 
1,2,-•• ,N) mostly equals zero except the frequency components really con¬ 
tained in s{t). 

We give the constraint conditions which should be satisfied and the data 
we need directly. Denoting p = fn/fs and M = {N — l)/2, we only need to 
sample at rate fs = 1/ At to obtain M samples si,S 2 , - ■ ■ , sm at time points 
£, provided that the parameters satisfy 

— TV is a prime number and N = Az + 3(z G Z"*"); 

— p is an arbitrary positive integer coprime to N; 

— C = {I ■ At : I = rri^ mod fV, m = 1, 2, • • • , M}. 

For instance, if V = 11, M = 5, we start to sample at the time t = 0, then we 
should get the samples at the time points lAt,3At,4:At,5At,9At. Noticing 
that p = /n /fs is simply required to be coprime to N, the sampling rate fs 
may be much lower than the Nyquist rate. Generally speaking, N does not 
always meet the first condition, we can select another satisfactory N' > N 
and a proper sampling rate. 

Then the detection method is explained as follows. Because 

exp (^j2TrfnlAtj = exp ^j27r '^"^^ = exp ^ 



(14) 


(15) 


Then 


y = Ax, 


(16) 
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where 

/exp(^ 2 £:^) exp(i 2 £:^) 

exp(i2£:^) exp(l2£:^) 

^exp(^^^:^)exp(i2£:ff^)... exp(i2£J^)y 

In (fTHll . y and A are known. If x can be solved from (ITT)ll . Cn and cor¬ 
responding to fn will be known. Obviously, x is a A-sparse vector. This is 
the same type of problem with ©• The coherence and RIP of A have been 
verified in Sect. [5J so the problem (1161) can be solved by common sparse re¬ 
covery algorithms. That is to say, after sampling in accordance with the above 
requirements, the frequencies and amplitudes of the harmonic components can 
be calculated. 


j27T'pl‘^'N ^ \ 

) 


exp( ^ 

exp(i22L£2=.Ar 


N 


(17) 


3.3 Practical Implementation 


In practical applications, real numbers are handled instead of complex num¬ 
bers. We give the implementation of our method in real number domain. For 
the convenience of the following discussion, we assume that the last element 
of X is 0, so the last column of A can be abandoned. Then (USD changes into 


y = 



(18) 


where B consists of the first M columns of A and B* is its conjugate, xi 
consists of the first M elements of x and X 2 consists of the second M elements 
of X adversely. By separating the real part and imaginary part, dT51) turns into 


Yr jyi — B7.X17. BjXij B^X2r B^X2j 

-(- j (B.^X]^^ -(- B,.X]^2 B^X27' -(- B,-X2j), 


(19) 


that is, 


= B'j 


where 





/ 


/ Xir\ 
Xli 


y^2r 

\^2^j 


,B' 


/ Br -Bj Br Bj \ 

Bj Br -Bi Br ) ' 


( 20 ) 


( 21 ) 


The subscripts r and i denote the real part and imaginary part of a matrix 
respectively. The compound matrix B' in (l20ll satisfies the same [k, ^)-RIP 
with A in real number domain. However, in practical measurement, such as in 
power systems, we couldn’t get the imaginary part y^. We solve this problem 
by assuming that X 2 of the actual signal is 0. That is to say, the highest 
harmonic frequency which can be detected is reduced by half. Comparing with 
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the theoretical analysis in complex number domain, we should choose twice 
the highest frequency as upper limit in practice. 

In the process of recovering x', we set X 2 r = and X 2 i = —Xii, namely 



( 22 ) 


V-xi*/ 


By using CS recovery algorithms, we obtain Xi^ and Xi^. Then the ampli¬ 
tude Cn and phase angle with respect to the discrete frequency /n(n = 
1, 2, • • • , M) can be determined. 

4 Simulation Results 

In this section, we shall create satisfactory signals and estimate their spec- 
trums by the proposed method. The signals are composed of several cosine 
waves. The frequency components are assumed to distribute randomly in 
{lHz,2Hz, - ■ ■ ,50Hz}. The amplitudes and phase angles of the frequency 
components are set to be uniformly distributed in [0.1,1] and [0, 27r). Accord¬ 
ing to our method, we set N = 103 and M = 51. The sampling frequency fs 
is set to be /at/IOO = 1.03Hz. The number of nonzero frequency components 
is known as the sparsity k. 

Firstly the situation without noise is considered. We use OMP as before 
to recover the signals and count the probability of success for varying k. The 
probability of success for every k is the result of 10000 trials. In Fig. [31 we 
can see the success probability is 1 when /c < 10, which is much smaller than 
the critical value shown in Fig. |3J This is because the actual sparsity of x' in 
solving process is no longer k. In fact, the locations of nonzero elements of x' 
have strong correlation because x' consists of a complex vector’s real part and 
imaginary part. If the recovery algorithm is modified to take advantage of this 
structural sparsity, better performance will be achieved. 

Then we estimate the spectrum of a frequency-sparse signal in white Gaus¬ 
sian noise with SNR=30dB. The signal has 10 frequency components and the 
corresponding amplitudes are 0.1,0.2,-- - ,0.9,1 respectively. The estimated 
spectrum is compared with original signal as shown in Fig. |4l The frequency 
components are all found with little error. 

At last we reconstruct the signal according to the amplitudes and phase 
angles we solve, even though it is not our aim. The reconstructed signal is 
plotted to compare with the original signal in Fig. |5l The reconstructed signal 
and the original signal are both sampled at lOOHz in time domain. The relative 
accumulated error of 50 points is about 3.85%. This indicates that the samples 
at the Nyquist rate can be obtained approximatively from under-sampled data. 

In above simulations, the sampling frequency has no influence on the re¬ 
sults, provided that p = Jn/Js is coprime to N. The errors decrease as the 
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Fig. 3 Success probability of estimation depending on sparsity in noiseless condition. 



Fig. 4 Comparison between original signal spectrum and estimated spectrum. 


sparsity k decreases and the SNR increases, which is the same as standard CS 
model. 


5 Conclusion 

In this paper, we construct a class of deterministic sensing matrices and verify 
their property as sensing matrices. Based on the special structure of the matri¬ 
ces, a deterministic under-sampled method of harmonic detection is proposed. 
Very few samples are required to estimate the frequencies and amplitudes of 
the harmonic components. Through theoretical analysis and numerical exper¬ 
iments, the feasibility and robustness of the method have been verified. 















































Title Suppressed Due to Excessive Length 


11 



Fig. 5 Comparison between original signal and restoring signal in time domain. 
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